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Abstract 

Background: Influenza is a well known and common human respiratory infection, causing significant morbidity and 
mortality every year. Despite Influenza variability, fast and reliable outbreak detection is required for health resource 
planning. Clinical health records, as published by the Diagnosticat database in Catalonia, host useful data for 
probabilistic detection of influenza outbreaks. 

Methods: This paper proposes a statistical method to detect influenza epidemic activity. Non-epidemic incidence 
rates are modeled against the exponential distribution, and the maximum likelihood estimate for the decaying factor 
k is calculated. The sequential detection algorithm updates the parameter as new data becomes available. Binary 
epidemic detection of weekly incidence rates is assessed by Kolmogorov-Smirnov test on the absolute difference 
between the empirical and the cumulative density function of the estimated exponential distribution with 
significance level 0 < a < 1. 

Results: The main advantage with respect to other approaches is the adoption of a statistically meaningful test, 
which provides an indicator of epidemic activity with an associated probability. The detection algorithm was initiated 
with parameter ko = 3.861 7 estimated from the training sequence (corresponding to non-epidemic incidence rates 
of the 2008-2009 influenza season) and sequentially updated. Kolmogorov-Smirnov test detected the following weeks 
as epidemic for each influenza season: 50 - 1 0 (2008-2009 season), 38 - 50 (2009-201 0 season), weeks 50-9 
(201 0-201 1 season) and weeks 3 to 1 2 for the current 201 1 -201 2 season. 

Conclusions: Real medical data was used to assess the validity of the approach, as well as to construct a realistic 
statistical model of weekly influenza incidence rates in non-epidemic periods. For the tested data, the results 
confirmed the ability of the algorithm to detect the start and the end of epidemic periods. In general, the proposed test 
could be applied to other data sets to quickly detect influenza outbreaks. The sequential structure of the test makes it 
suitable for implementation in many platforms at a low computational cost without requiring to store large data sets. 
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Background 

Influenza is a well known and common human respira- 
tory infection. It is responsible of significant morbidity 
and mortality every year. The World Health Organization 
(WHO) estimates that annual epidemics result in about 3 
to 5 million cases of severe illness and about 250,000 to 
500,000 casualties worldwide [1]. Influenza has a seasonal 
pattern and in the Northern hemisphere epidemics occurs 
between November and March. 
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Sentinel networks covering less than 2% of the pop- 
ulation have been the traditional surveillance system. 
More recently, electronic health records are widely imple- 
mented in some regions making available a significant 
amount of health related data. In Catalonia, primary care 
doctors have been routinely registering their activity in 
eCAP (an electronic health recording system) since 2006. 
This accounts for over 3,500 physicians collecting data of 
nearly 6 million people (80% of the population) [2,3]. The 
data produced is a key source of information that could 
also be used for surveillance of certain diseases. In 2011 
SISAP (the Catalan acronym for Information Systems for 
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Primary Care Services) developed a process of extract- 
ing health conditions data from eCAP in order to provide 
information about some diseases. Nowadays, SISAP pub- 
lishes weekly information about all cases of those infec- 
tious diseases, available on Diagnosticat (http://4.sisap. 
cat/diagnosticat). 

Among the large variety of tracked diseases, we focus 
on influenza. Influenza data on Diagnosticat has shown 
its validity as compared to the sentinel network data. Its 
main advantage is that its data is available faster [3], and 
is thus theoretically better suited for timely activation of 
emergency responses. More precisely, we are interested 
in the outbreak detection of such disease by means of 
a statistical surveillance system [4]. The objective is to 
provide an indicator as soon as data gives enough evi- 
dence to assess the event. Moreover, we seek to apply 
the developed algorithms to the Diagnosticats database. 
Other influenza surveillance systems reported in the lit- 
erature are briefly commented in the following. In [5], 
the authors presented an online tool to analyze epidemi- 
ological data. They used non-epidemic training data to fit 
a time-series model for the periodic baseline level, thus 
following a parametric approach for their purposes. [6] 
proposed a method to detect deviations from the baseline 
that was based on fitting historical data to a time-series 
regression model. Time-series analysis were also used in 
the same context in [7] . A weighted moving average algo- 
rithm was considered in [8] to monitor the time-series. In 
contrast, in the present work we use non-epidemic data 
to infer the parameters of the probability density function 
of non-epidemic periods. Therefore, we make no assump- 
tions on the time evolution of the disease and instead we 
resort to a probabilistic characterization. The work in [9] 
suggests that non-parametric methods are best suited to 
influenza modeling, due to its large variability between 
years. The authors of [9] focused on other features of the 
influenza cycle rather than disease outbreak, such as the 
peak after which the influenza incidence starts to decline, 
and applied the method to a case study in Sweden. More 
involved methods, combining both approaches, are pos- 
sible. In [10] a switching Markov model was used with 
an autoregressive process to model epidemic data and 
a white Gaussian process for non-epidemic modeling. 
This method was implemented in a web-based applica- 
tion developed by the same authors and presented in [11]. 
Other references point to imaginative, yet reliable meth- 
ods to detect possible epidemic outbreaks. For instance, in 
[12], Google search queries were used to track influenza- 
like illness in a population. In a similar fashion, Twitter 
was used in [13]. 

This paper presents a statistical surveillance system 
which provides an automated detection of influenza epi- 
demics. Health resource planning is an application which 
could benefit from this tool. The proposed method is 



able to operate on-line a and it is based on the statisti- 
cal characterization of non-epidemic influenza incidence 
rates. A major advantage of this approach is its statistical 
meaningfulness, and thus detection is not only a binary 
result but the confidence in its outcome can be assessed 
in terms of probabilities resorting to hypothesis testing 
theory. Non-epidemic data is modeled with an exponen- 
tial distribution, in the vein of [14]. The unique parameter 
of the distribution is the decaying factor, that is initially 
estimated by a training data set and sequentially updated 
as new observations are recorded. Such statistical char- 
acterization is used to design a detector based on the 
one-sample Kolmogorov-Smirnov test. The method was 
applied to data in the Diagnosticat database to successfully 
detect influenza activity. 

The remainder of the paper is as follows. Diagnosticat 
is introduced in the following section Then we provide 
insights on the statistical distribution of influenza inci- 
dence rates, as well as how the relevant parameters can be 
estimated from the observations. The general statistical 
detector is proposed, and the results for the Catalan case 
study are presented. 

Methods 

Diagnosticat: an open epidemiological database 

Diagnosticat is an open-access database which contains 
reports of many diseases occurring in Catalonia, such 
as influenza, papilloma or chickenpox. The informa- 
tion available in the Diagnosticat database includes all 
clinical influenza diagnoses codes (ICD-10 code) and 
is obtained weekly from eCAP through an automated 
process. The website is timely updated a few minutes 
after every finished epidemiological week (EW). After 
the extraction, a computer algorithm automatically cre- 
ates the different tables with the information that is used 
in the website. No identifiable or personal information 
on patients is used, maintained or transferred through 
this system. 

Currently, the Diagnosticats database is composed of 
data from 4 influenza seasons. Information is available 
starting on 2008 and is updated weekly since then. Data 
is presented as incidence rates per 10 5 population, a unit 
that allows comparison of diagnoses over different terri- 
tories independently of the number of inhabitants. In this 
work, the EW is a group of seven days that begins on a 
Sunday and ends on a Saturday. 

Although communicable diseases are in general yearly 
represented, influenza is represented by seasons due to the 
characteristics of the influenza epidemic. The epidemic 
usually starts at the end of the year and ends mid of the 
following year, peaking in December and January. For this 
reason and to be consistent with the influenza epidemic, 
Diagnosticat uses graphics by seasons beginning the EW 
23 and ending at EW 22 the following year. 
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Statistical data analysis 

In this section we aim at obtaining a statistical model for 
the recorded data of influenza cases. Particularly, we find 
out that an exponential distribution might be an appropri- 
ate way of modeling non-epidemic data. Let us consider 
that the rate of influenza cases in non-epidemic periods is 
a random variable X and that a database contains weekly 
realizations of such variable. We define xt as the set con- 
taining chronologically-ordered observations up to time t, 
X t = {xi,...,Xt) e Let us consider that, under the 
hypothesis that influenza is not active, X is exponentially 
distributed [14], with probability density function (PDF) 
given by 



fx(x;X) = Xe 



-Xx 



x>0 



(1) 



and 0 otherwise. Mean and variance are expressed in 
terms of the parameter X > 0 as K{X} = 1/X and 
Var{X} = 1/X 2 , respectively [15]. Recall that only a sub- 
set of xt follows (1), that is those observations taken in 
non-epidemic weeks. We define this subset as & c Xt and 
its elements {x\, . . . ,XL t ] are those from xt which were 
detected (using the herein proposed method for instance) 
as in non-epidemic periods. We use L t to denote the total 
number of non-epidemic weeks up to t. 

Such statistical characterization means that, if not in an 
epidemic scenario, the rate of cases would be most likely 
close to zero and decreasing according to the exponential 
factor A. The factor X is unknown and thus has to be esti- 
mated for a complete characterization of (1). In the sequel, 
we derive the Maximum Likelihood (ML) estimator of X 
given the random set Assuming that observations are 
independent, the likelihood function of X given & is 



C(X) = X Lt e~ x ^=i~ Xl , 



(2) 



whose optimization is equivalent to maximization of the 
log-likelihood. Derivative of the latter yields to 



d\nC(X) 
dX 



U 
X 



u 

i=i 



(3) 



and equating to zero we obtain the ML estimate of the 
exponential factor as 



1 i=i / 



(4) 



which is unbiased and asymptotically consistent with vari- 
ance Var{A^} ~ Var{X} 3 /Z,£. Notice that X t is computed 
using the L t observations in the non-epidemic period up 
to the current week t. 

We study the goodness-of-fit of the statistical model in 
our case study, i.e. the database in Diagnosticat. Figure 1 
shows a comparison of the histogram of non-epidemic 
influenza incidence rates and the assumed exponential 
distribution, with X = 4.6899 estimated as in (4). In order 



to provide more observations to evaluate the goodness of 
the exponential fitting, in this case we used all records in 
the database and omitted those observations that where 
considered a posteriori in the influenza epidemic phase. 
To identify these data we have selected as the epidemic 
phase the same weeks as those established by the sen- 
tinel surveillance system in Catalonia [16], refer to row 
Sentinel network in Table "Start and end of influenza epi- 
demic in EW per season as detected by the method" to 
identify the periods. However, in an automated surveil- 
lance system, X should be estimated sequentially with an 
initial guess computed using records from the first sea- 
son, which could be used as a training sequence. We also 
run a/ 2 goodness-of-fit test to validate the exponential 
assumption, obtaining a p- value of 0.8451. 

Sequential influenza detector 

The objective is to build an autonomous detection algo- 
rithm that is able to determine whether influenza is active 
or not based on the records of influenza cases. In the case 
of Diagnosticat, these observations are weekly received, 
and thus a week-by-week detection is provided by the 
method. As mentioned in [17], such an automated surveil- 
lance system needs to fulfill two objectives: timeliness, to 
ensure that intervention protocols to the epidemic are in 
time; and sensitivity, to avoid time/effort wasting due to 
precipitated measures. 

The detector is based on the idea that we have a sta- 
tistical characterization of the influenza incidence rates 
when not in the epidemic phase. As discussed earlier, 
these observations follow an exponential distribution with 
parameter X estimated as in (4) in batch processing or as 
in (8) for the sequential treatment of data. 

For a given instant t, the null hypothesis Tio is that 
the observation Xt comes from a distribution as speci- 
fied in (1) with Xt-i, being the ML estimation at time 
t — 1 of the exponential factor. If the hypothesis Ho is not 
accepted, it means that influenza is detected at the £-th 
EW. Equivalently, we say that 



% : x t ~fx(x;X = X t -i) . 



(5) 



We notice that the test has to be conducted using a sin- 
gle observation, x t , from the random variable and that the 
baseline distribution is continuous and completely speci- 
fied. Although other alternatives might apply, in this type 
of decision problems we could obtain enhanced perfor- 
mance by resorting on Empirical Density Function (EDF) 
statistics to assess Ho, as claimed in [18]. Since we have 
only one observation available, the EDF is directly Fx (x) = 
l{x t <x}> where 1{A} is the indicator function of an event A, 
equal to 1 if A is true and 0 otherwise. Here we use the 
classical Kolmogorov-Smirnov (KS) nonparametric test 
[19], which is based on the maximum absolute difference 
between the EDF and the Cumulative Density Function 
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100 



(CDF) of the reference distribution. The latter is known to 
be Fx(x; X) = 1 — e~ Xx for the exponential distribution. 
In general, the one-sample KS statistic is denned as 



D (Px(x),Fx(x)^ = sup \P x (x) - F x (x)\ 



(6) 



although in our case the statistic D can be simplified 
accounting that only one observation x t is used to con- 
struct the EDF, and thus the EDF can only take values 0 or 
1: 



D (Px(x),F x (x)J = max {l 
■F x (x)\ 



-Xx 



-Xx\ 



(7) 



0- 1 + e 



-Xx I 



where we used that |0 
e -Xx w h en p x (x) = 0. 

With this detector, hypothesis Ho is rejected at a signif- 
icance level 0 < a < 1 if D > y a , where the threshold 
satisfies a = P {D > y a }. Relevant threshold values can be 
consulted from tables [20,21] if Fx(x) is completely spec- 
ified, and from [22] if its parameters should be estimated 
from data. Notice that the latter is the case here, since k is 
estimated. 

The result of the hypothesis test is binary. Another use- 
ful way of reporting the result of the test are ^-values. 
Given d, a realization of the random variable D, we recall 
that the j?-value is defined as the probability of obtain- 
ing another realization of the test at least as extreme 
as d conditional on Ho being true, mathematically p = 
P {D > d\ Ho} (or Type I error). Therefore, not only Ho is 
rejected if the p-vdlue of the test is less than a, but p is 
also quantifying the confidence in such decision. Large p- 
values indicate that xt is likely to be generated from Fx(x), 
and viceversa. 



To complete the sequential detection algorithm, we 
need to estimate the exponential factor k of the baseline 
distribution. Here we are interested in sequential treat- 
ment of epidemic data, and therefore we would like to 
use an alternative expression to the batch approach in (4) 
that allows us to compute it recursively instead of pro- 
cessing the entire set of L t non-epidemic weeks each time 
a new observation is recorded. The propose procedure is 
applied after the current observation x t is processed by 
the KS-based detector. If based on xt the method rejects 
Ho> then we keep the same estimation as previously, k t = 
kt-i. Otherwise, if x t shows evidences that at the £-th EW 
influenza was active, then we count another week in L t 
and use the following recursive expression 



a: 



(L t -l)+x t \ 



(8) 



which is algebraically equivalent to (4). Using (8) instead 
of (4) has the advantage that new data is processed upon 
arrival, and thus there is no need to store the complete 
dataset nor reprocess all data for each new measurement. 
Initialization of the exponential factor Xq could be com- 
puted for instance based on expert, a priori knowledge. 
In our implementation of the method, we used the non- 
epidemic weeks in the first influenza season as a training 
sequence to estimate the value of ko. L t counts the num- 
ber of times the null hypothesis was assessed valid up to 
t, and thus it is also the amount of observations used to 
estimate k. 

The pseudo-code of the sequential influenza detector 
can be consulted in the Algorithm. Recursive exponen- 
tial fitting has also been included, note that k t is updated 
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only when Ho is accepted. For the sake of clarity, Figure 2 
depicts a block diagram representation of the method. To 
sum up, when a new observation x t is recorded, it is tested 
against the exponential distribution and a new value of 
the statistic D is computed which is used to accept or 
reject Ho- The parameter of the exponential distribution 
at the null hypothesis is computed based on all previous 
observations which were not in epidemic phase, and thus 
we write k t -i. This estimate is only update according to 
(8) if Ho is accepted based on the current observation, 
otherwise the previous estimate is kept. 



Algorithm Sequential detection of influenza epidemics 

1: Initialization t = 1, An, and Lo = 1 

2: At time t new data becomes available, x t 

3: Compute statistic D = max 1 1 - e -^t-m 

4: if D > y a then 

5: Reject the null hypothesis =>► Flu detected 

at the Nth EW. 
6: Keep k t = k t -i and L t = L t _\ 

7: else if D < y a then 

8: Accept the null hypothesis => Flu not 

detected at the t-th EW. 
9: Set L t = L t -\ + 1, the number of 

observations used to calculate k t 
10: Update k~ l = ± (k~\(L t - 1) + x t ) 

11: end if 

12: t=t+l and go to step 2 



Results and discussion 

We used the open database described earlier to test the 
detector. We also validated in this experiment the pro- 
posed ML data fitting that characterizes non-epidemic 
cases as exponentially distributed. At this time of writ- 
ing, we had available data from the 2008-2009 to the 
2011-2012 seasons. 

The parameter k was estimated sequentially as in (8). 
Initially, this value was calculated using data from the first 
influenza season (i.e., 2008-2009), as a training sequence 
for ko, with those observations identified by the Cata- 
lan sentinel network as being in the epidemic phase 
[16] removed from that years set x to form the sub- 
set f . In particular, incidence rates per 10 5 population 
above 20 were not considered. If not otherwise stated, we 
used 20/ 10 5 to train the method, although later in this 
section we provide a sensitivity analysis with respect to 
this threshold. The resulting exponential factor was ko = 
3.8617 for the first year, which was updated as described in 
the algorithm when new data from EW 23 of the influenza 
season 2009-2010 was available. 

We tested the observation corresponding to the £-th EW 
according to (7) with a significance level a = 0.05. The 
results of the test (red crosses) and the registered influenza 
incidence rates (blue solid line) used to infer the decision 
are plotted in Figure 3. Red crosses denote the EW where 
Ho was rejected, i.e. influenza activity was detected. The 
test detected the following weeks as epidemic for each 
influenza season: 50 - 10 (2008-2009 season), 38 - 50 
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Figure 3 Weekly influenza-detection results for 4 influenza seasons along with recorded data, y-axis represents the influenza cases per 
10 5 inhabitants (solid blue) and the detected influenza outbreaks output (red crosses). 



(2009-2010 season), weeks 50 - 9 (2010-2011 season) 
and weeks 3 to 12 for the current 2011-2012 season. The 
results obtained were in accordance to the assessment 
provided by the official system of influenza surveillance 
in Catalonia (PIDIRAC) [16,23-25]. PIDIRAC establishes 
the epidemic periods after the analysis of data from sen- 
tinel network and virus activity (laboratory surveillance 
data). In Table 1 we summarize the results in terms of 
start and end of influenza periods per season. Comparing 
the obtained results to those given by the official source 
we can calculate for the method an accuracy of 90%, a 
sensitivity of 100% and a specificity of 87.6% for the com- 
plete four seasons. Recall that we used a threshold on the 
first season of 20/ 10 5 to estimate Xq based on [16]. The 



selection of this value needs to be done depending on the 
application. For the sake of completeness, the table also 
shows a sensitivity analysis to such threshold. The con- 
clusion is that if this parameter is set within reasonable 
intervals (as discussed earlier, based on prior knowledge), 
the results of the detector are not drastically varying. 
Indeed, the most sensitive season is the first, which is 
used for Xq calculation. Detection results for other sea- 
sons remain unaltered. Increasing the threshold makes the 
method less sensitive to epidemic data, and viceversa. 

It is important to notice that in the influenza sea- 
son 2009-2010, the A(H1N1) Influenza virus pandemic 
occurred during autumn in Catalonia, fact that caused 
a different temporal pattern of influenza epidemics [26]. 



Table 1 Start and end of influenza epidemic in EW per season as detected by the method 





2008-2009 




2009-2010 




2010-2011 




2011-2012 




Threshold 


Start 


End 


Start 


End 


Start 


End 


Start 


End 


1 0/1 0 5 


50 


11 


38 


50 


50 


9 


3 


12 


20/10 5 


50 


10 


38 


50 


50 


9 


3 


12 


30/1 0 5 


51 


8 


38 


50 


50 


9 


3 


12 


40/10 5 


52 


7 


38 


50 


50 


9 


3 


12 


50/1 0 5 


52 


7 


38 


50 


50 


9 


3 


12 


Sentinel network 


51 


8 


41 


51 


51 


9 


4 


12 



Sensitivity to a number of thresholds on the first season to estimate Xq is reported. Also, the official reports based on the sentinel network are included. 
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2012 



Recall that similar differences were reported worldwide 
[27,28]. That epidemic had higher incidence rates than 
seasonal influenza and took place some weeks before. The 
proposed method was also able to cope with this non- 
seasonal event as its departure from Ho was similar to that 
of seasonal influenza. 

Figure 4 shows the evolution over time of X t , 
where the initial training can be observed as the 



constant value A, = 3.8617. It becomes clear from this 
figure that the sequential updating is only performed in 
non-epidemic periods, and constant otherwise. For the 
sake of completeness, we compared the results of the 
sequential method to the estimates obtained after pro- 
cessing all data using the batch method in (4). Notice 
that both methods provide the same estimate after all 
data is processed, with the advantage that the sequential 
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Figure 5 Sequence of p-values obtained by the detector. 
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procedure does not require storage of all observations 
before processing them. 

In order to provide more insight in the behavior of the 
detector, Figure 5 shows the obtained p- values. As com- 
mented earlier, this value quantifies how likely the null 
hypothesis is true. We can observe that low values were 
obtained when Ho was rejected and that it grew in non- 
epidemic periods. In the latter, ^-values were not very 
large (ideally close to one) due to very small sample size, 
recall that we were constrained to use a single observation 
per EW to infer the decision. It is remarkable that, despite 
the small sample size, results were in accordance to what 
the official system based on the sentinel network recorded. 

Indeed, the proposed method can be used in general to 
the detection of other infectious diseases whose statisti- 
cal characterization is available. If non-epidemic periods 
could be characterized by an exponential distribution, 
the usage of the method is straightforward. Otherwise, if 
another distribution better fits the data, slight modifica- 
tion of the method should be performed to compare it 
with the EDF and to update the required parameters of 
the PDF. A limitation of the method is related to the data 
gathering method, which has to be continuously recording 
observations. Since the method is based on non-epidemic 
data to estimate the distributional parameters and to 
assess whether the disease is active or not, some systems 
like sentinels cannot straightforwardly benefit from this 
tool. Recall, for instance, that sentinels are likely to stop 
recording data in typically inactive periods of the disease, 
these inactive periods have an impact in the quality of the 
estimated distribution of non-epidemic data. 

Conclusions 

In this paper we proposed an automated method to detect 
influenza outbreaks from periodically recorded incidence 
rates. In contrast to setting yearly predefined thresholds 
to determine influenza outbreaks by data inspection, we 
presented a detector based on the statistical properties of 
non-epidemic data. The method can be useful to com- 
plement traditional surveillance methods. The algorithm 
provides a binary signal indicating epidemic activity as 
well as a quantitative measure (i.e., the p-vdlue) of its con- 
fidence. The method can be executed sequentially and it 
is self-adjusted (after some initialization of the exponen- 
tial factor). We presented results with real medical data 
from Catalonia, showing the performance of the detec- 
tor even in non-seasonal scenarios, as well as the validity 
of the exponential modeling for non-epidemic influenza 
incidence rates. Finally, the objectives required for such an 
automated system are seen to be fulfilled: timeliness and 
sensitivity. 

Timeliness is generally defined as the difference between 
the time an event occurs and the time the reference stan- 
dard for that event occurs. Diagnosticat data is available 



the instant the epidemiologic week has endend, and thus 
advances in four days the publication of data over the sen- 
tinel network based system in Catalonia. Additionally, the 
proposed sequential detection method signals the event 1 
week earlier on average for the tested data (excluding the 
AHlNl.pdm.2009 season). 

Endnote 

a Note that in this work we use on-line to denote that 
the algorithm operates sequentially, that is, as new data 
becomes available. 
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